Thermodynamics of Heisenberg ferromagnets with arbitrary spin in a magnetic field 
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The thermodynamic properties (magnetization, magnetic susceptibility, transverse and longitu- 
dinal correlation lengths, specific heat) of one- and two-dimensional ferromagnets with arbitrary 
spin S in a magnetic field are investigated by a second-order Green-function theory. In addition, 
quantum Monte Carlo simulations for S = 1/2 and S = 1 are performed using the stochastic series 
expansion method. A good agreement between the results of both approaches is found. The field 
dependence of the position of the maximum in the temperature dependence of the susceptibility 
fits well to a power law at low fields and to a linear increase at high fields. The maximum height 
decreases according to a power law in the whole field region. The longitudinal correlation length 
may show an anomalous temperature dependence: a minimum followed by a maximum with in- 
creasing temperature. Considering the specific heat in one dimension and at low magnetic fields, 
two maxima in its temperature dependence for both the 5 = 1/2 and S — 1 ferromagnets are found. 
For S > 1 only one maximum occurs, as in the two-dimensional ferromagnets. Relating the theory 
to experiments on the S = 1/2 quasi-one-dimensional copper salt TMCuC [(CHa^NCuChj], a fit to 
the magnetization as a function of the magnetic field yields the value of the exchange energy which 
is used to make predictions for the occurrence of two maxima in the temperature dependence of the 
specific heat. 

PACS numbers: 75.10.Jm, 75.40.Cx 



I. INTRODUCTION 

The study of low-dimensional quantum spin systems^ 
is of growing interest and is motivated by the progress 
in the synthesis of new materials, where ferromagnetic 
compounds attract increasing attention. For exam- 
ple, besides the spin S = 1/2 quasi-one-dimensional 
(ID) ferromagnetic systems, such as the copper salt 
TMCuC^ the organic magnets p-NPNNi£ and [3- 
BBDTA-GaBr4j£ and the CuCl2-sulfoxide complexes^ 
recently the S = 1/2 quasi-2D ferromagnet Cs2AgF4, 
which has a structure similar to the high-T c parent 
compound LagCuCU, was studied 8 and found to be 
magnetically reminescent of K2CUF412 In ferromagnetic 
systems with S ^ 1 mainly the effects of single-ion spin 
anisotropies were investigated, such as in the quasi-lD 
5 = 1 easy-plane ferromagnet CsNiF^ and in 2D 
easy-axis Heisenberg models in a magnetic fiel d 11 ' 12 ! 13 
for describing the spin reorientation transition in thin 
ferromagnetic films (see Ref. [13 and references therein) . 

The 2D anisotropic S ^ 1 Heisenberg ferromagnets 
in a magnetic field were investigated by Green-function 
methods ) 11 ' 13 ' 15 where the exchange term was treated 
in the random-phase approximation (RPA)^ and by 
quantum Monte Carlo (QMC) simulations. 12 In a pre- 
vious paper— we have developed a second-order Green- 
function theory of ID and 2D 5 = 1/2 ferromagnets in 
a magnetic field which goes one step beyond the RPA 
and provides a rather good description of magnetic short- 
range order (SRO) and of the thermodynamics. This can 
be seen from the comparison with the exact calculations 
by the Bethe-ansatz method for the quantum transfer 



matrix in the ID model and with the exact diagonaliza- 
tions on finite lattices. In particular, for the 5=1/2 
ferromagnetic chain two maxima in the temperature de- 
pendence of the specific heat at very low magnetic fields 
were found. On the contrary, the RPA was shown to 
fail in describing the SRO, reflected, e.g., in the specific 
heat, whereas the magnetization and the magnetic sus- 
ceptibility are quite well reproduced. Recently, a similar 
Green-function approach for 5=1/2 ferromagnets was 
presented^ which improves the theory of Ref. [13 con- 
cerning the agreement with exact methods. The results 
obtained for 5=1/2 are stimulating to investigate fer- 
romagnets with 5 > 1/2 in a magnetic field, for which 
a second-order Green-function theory of SRO is not yet 
developed. Second-order Green-function approaches for 
ferromagnets with arbitrary spin exist in the case of zero 
magnetic field onl y 19 ' 20 where in Ref. [2^ ferromagnetic 
chains with an easy-axis single-ion anisotropy were stud- 
ied. 

In this paper we extend both our previous theory for 
5 = 1/2 (Ref. [13) to arbitrary spins and the theory of 
Ref. [l9| for zero field to arbitrary fields. We start from 
the ferromagnetic Heisenberg model with arbitrary spin 
5, 



H = -jJ2SiSj-hJ2Si 



(i) 



denote nearest-neighbor (NN) bonds along a chain 
or on a square lattice; throughout we set J = 1] with 
Sf = 5(5 + 1). We calculate thermodynamic prop- 
erties (magnetization, magnetic susceptibility, correla- 
tion length, specific heat) at arbitrary temperatures and 
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fields. For comparison, we perform QMC simulations 
of the S — 1/2 and 5=1 models on a chain up to 
N = L = 1024 sites and on a square lattice up to 
N = L x L = 64 x 64. 

The rest of the paper is organized as follows. In Sec.HTl 
the second-order Green-function theory for model (JTJ) is 
developed, where the extensions of previous second-order 
Green-function approache a 17 ' 18 i 19 to arbitrary spins and 
fields imply novel technical aspects. Moreover, consider- 
ing the case 5=1/2 the theory is extended as compared 
with Refs. [l?] and [H by the introduction of two addi- 
tional vertex parameters and, correspondingly, by taking 
into consideration two additional conditions for their de- 
termination. This extension is shown to have qualitative 
effects on the temperature dependence of the longitu- 
dinal correlation length (see Sec. IIVI B). In Sec. IIIII the 
employed QMC method is briefly described. In Sec. Hvl 
the thermodynamic properties of the ID and 2D ferro- 
magnets are investigated as functions of temperature and 
field, also in comparison with RPA, and are related to 
experiments. Particular attention is paid to the calcula- 
tion of the transverse and longitudinal correlation lengths 
which were not considered in Refs. Qj] and [TH Finally, a 
summary of our work is given in Sec. |VJ 



II. SECOND-ORDER GREEN-FUNCTION 
THEORY 

To determine the transverse and longitudinal spin cor- 
relation functions and the thermodynamic quantities, we 
employ the equation of motion method for two-time re- 
tarded commutator Green functions^ First we calcu- 
late the transverse spin correlation functions. Because 
we treat arbitrary spins in nonzero magnetic fields, so 
that we have (5 Z ) ^ 0, we consider the Green functions 

((5+; S_ q )) u introduced by Tyablikov within the first- 
order theory, i.e. the RPA (see Appendix), where S_l 
is the Fourier transform of = (5f) n 5 i ~ with n = 

0, 1, 25 — 1, and the Green functions ((«5+; S_g)) u 
which we calculate for the first time in the second-order 
theory. The equations of motion read 

c«5+; S™-)) u = M(")+" + «i5+; S^~)) u , (2) 

u>((iS+; SM-)) U = + «-5+; S™~)) u . (3) 

The moments M^ + ~ = ([S+,S^~]) and M q n)+ ~ = 
([iS q , S_ q }) are given by the exact expressions 

M^+- = 2((5T +1 } + (1 - «»,o) E it) C" 1 )** 

k=l ^ ' 

{S(S + l)((S z ) n ~ k ) + ((S z ) n - k+1 ) - ((S z ) n - k+2 )}, (4) 



Mf - = zil J q ){2C^ )zz + C^- + 

+ (l-5„,o)EQ(-l) fe [5(5 + l) 

x (S k , n (S z ) + (1 - 8 k , n )C§- k - l)zz ) 

+ C^ k)zz -C^ k+1) "]} + hM^ + -, (5) 

where Cnm = = (5q ^ ^r)> C n m = 

z/2 

Cr )ZZ = ((S§) n+1 S Z R ), R = ne x +me y , lq = § ^ cos ft, 

i=i 

and z is the coordination number. Deriving Eqs. ([4]) and 
([5]) the operator identity 

Sf=SrS+ + S z + (S z ) 2 (6) 

has been used. In Eq. the second derivative — S q 
is approximated as indicated in Refs. 11?. 19.20. 21 l22ll23 . 
[241 . That means, in —5^ we decouple the products of 
operators along NN sequences as 

S+S+S- = a+- (SfSf)S+ + a+- (S+Sf)S+, (7) 

where the vertex parameters a\ and oi\~ arc attached 
to NN and further-distant correlation functions, respec- 
tively. The products of operators with two coinciding 
sites, appearing for 5 ^ 1, are decoupled a o 19 ' 20 

St$j Sj = (Sj 5^)5+ + A + (Si Sj )Sj , (8) 

where the vertex parameter is introduced. We ob- 
tain 

-S+ = [(cj+-) 2 -h 2 }S+ + 2hiS+ (9) 

with 

«~) 2 = f (1 - 7q){A + - + 2za+-C w (l 7 q)}, (10) 

A+- = 5(5+l) + <(5*) 2 ) 

+ 2{A+--(z + l)a+-}C 10 

+ 2a+-{(z-2)C 11 + C 20 }, (11) 

where C nrn — ^C n °m + + Cnm Z ■ In the special case 
5 = 1/2, in — 5j + products of spin operators with two 
coinciding sites do not appear which is equivalent to set- 
ting A H = 0. Finally, we get the Green functions 

«5+;5W->>.= E ( 12 ) 



i=l,2 q 
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where 

u q i t2 =h±u+-, (14) 

4% = \m^ + ~ ± ^(M(»)+- - MfW+") (15) 

with the moments given by Eqs. ((H) and |5]). The 
transverse dynamic spin susceptibility x+~(o;) = 
-((5+; £:,))„, is given by Eq. for n = 0. 

Because we consider nonzero magnetic fields within the 
second-order theory, the behavior of the Green functions 
(fT2")l with the poles fT4")) exhibits, for arbitrary spin, a 
peculiar aspect. Considering the static Green functions 
({S£;S_q ))w=0) in particular the static spin suscepti- 
bility Xq ~ = Xq~( w = 0), a divergency signaling a phase 
transition could appear if uj q 2 — 0, i.e., = h. Ac- 

cording to Eq. IjlOp the corresponding q values are given 
by 

1-7, = go = {Aza^CwY 1 x 

{[(A+-) 2 + 16 a +-C 10 /i 2 ] 1/2 - A+-}. (16) 

This equation may be fulfilled in region I of the h — T 
plane defined by h < ho(T), where ho(T) is determined 
by Eq. (fl6|) with go — 2 which is realized at the corner 
of the Brillouin zone with 7, = — 1. In region II, h > 
ho(T), we have h > for all q. For nonzero fields 

the Heisenberg ferromagnet described by Eq. ([T]) has no 
phase transition. This means, Xq ~ has to be finite at all 
q. We require this regularity to hold also for the static 
Green functions with n = 1, 25— 1. That is, in region 
I we require — with q given by Eq. (p~6|) . This 
results in the regularity conditions, 

MfW+- = z{2C[f zz + C$~ + + 

V - § ^ £ (fc) ("^"^(^ + Wk,n{S Z ) + 

{ -1 c \^f(ri— fe— l)zz\ /-^{n—k^zz /~<(n— fe+l)z«i-i 
(l-0fe, n jO 10 J + -°io J/So- 

(17) 

Note that Eq. (TT7)) for 5 = 1/2 agrees with the condition 
given in Ref. [l8| which is obtained from an analyticity ar- 
gument and is written as an expression for (S z ). In the 
limit T — > 00, the field ho separating the regions I and II 
may be easily obtained. For T — > 00 we have spin rota- 
tional symmetry so that <(5 Z ) 2 ) = \C$~ + . By Eq. @ 
with lim (S z ) = we get ((5 2 ) 2 ) = ±5(5 + 1) result- 

ing in A+- = f 5(5 + 1) and {uj+-f = §A+"(1 - 7g ). 
From = h and go = 2 we get lim ho(T) — 

2^zS(S + l)/3. Following Ref. [H, we assume the con- 
ditions (fl7|) to be valid also in region II. This guarantees 
the continuity of all quantities at the boundary ho(T). 
From the Green functions (| 1 2[) and (|13p the 

transverse correlators C% } ~ + = (l/iV^C^"-^^ 

q 



and C^ -4 " = (l/N)J2C { q n) ~ + e iqR with the struc- 
q 

ture factors C q n) - + = (S^~ 5+) and C< n) - + = 
(5^~«5q) are calculated by the spectral theorem, 

i=l,2 

<5<»)-+ = S^V^), 

i=l,2 

where n(w) = (e^ - l)" 1 and 0=1/T. 

Now we derive some useful sum rules. Using 
(5| n) ~5+) = ({Sf) n SrS+) obtained from Eq. © multi- 
plied by (5?)™ (n = 0, 1, 25 - 1) and Eq. |jlgj) we get 
the relation 

5(5 + 1)((5 2 )") - ((5 Z )" +1 ) - ((5 2 )"+ 2 > 

4EE4 )fl W' (19) 

q i=l,2 

S 

By the identity 1 | (5? — m) = one can express 

m=-s 

(Sf ) 2S+1 appearing in Eq. (flU for n = 25 — 1 in terms 
of lower powers of Sf (Refs. [lg andl25l). 

2S 

(S~) 25+1 = $>L 5) (5f) fe , (20) 

fe=0 

where the coefficients otjp are given in Ref. [25]. From 
the system of the 25 equations ([IT?]) we can determine 
the magnetization m = — 2/ib(S z ). 

Similarly, in the second-order theory higher-derivative 
sum rules may be derived which, for nonzero fields, pro- 
vide 25 additional equations for determining the ver- 
tex parameters and some longitudinal correlators (see 
below). Multiplying S$ n) ~ by iSf = ^ [S Z S+ - 

j(n.n.i) 

5/5?) + hS+ and using Eqs. ((T3J), JH}, dHJ) and © 
we obtain 

z{5(5 + 1)[<5„, <5 2 ) + (1 - t^cft-V"] 
— u 10 - o 10 - o 10 - o 10 J- 

= -^E ( " 1)<w i"^ BW - (21) 

q i=l,2 

The correlator c}o +1 ^* for n = 25— 1 may be expressed 
in terms of (S z ) and c[q )zz with n € 25-1 by Eq. ([20|. 

Equally, C 1 2 | S ' ) h can be written in terms of Cy^ h (n ^5 
25 - 1) by the identity^ 

2S-1 

Sr(S z r s = Sr £ 4 Sll) (5?) fe , (22) 
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where the coefficients S/f are given in Ref. Hfl The 

product (S Z ) 2S S~ appearing in c[q S ^ h can be deduced 
from Eq. (|22| by the commutation relations for spin op- 
erators. The sum rule (f2"T|) for n — also follows from 
the exact representation of the internal energy per site, 
u = (H)/N = -f (C[l ] ~ + + C§ )zz ) - h(S z ), in terms of 
((S q ; SZ q ))u which can be derived similarly as in Ref.l26l 
for 5* = 1/2, 



^[5(5 + 1)<5*) - c[ )zz - - h(S>) 



N ^ 



2^ 



(u-h)Im((S+;SZ q )) u n(u), ( 2 3) 



if the result (TT2|) for ((S q ; S_ q )) u (n = 0) is inserted into 
Eq. ((33}. 

To calculate the longitudinal spin correlation functions 
C^ zz from the Green function {(S q ; 5 2 q )) w = — x q 2 (w), 
where x q z (^) is the longitudinal dynamic spin suscepti- 
bility, we start from the equations of motion analogous 
to Eqs. |(5J| and ([3} and perform a second-order decou- 
pling which is equivalent to the projection method with 
the basis (S q ,iS q ) neglecting the self-energy, as indi- 
cated in our previous papers * 17 ! 20 In —Sf we adopt the 
decoupling a 17 ' 19 i 20 analogous to Eqs. {7} and ©, 



SI 5^5; — ol\ (Sj S[ )S i , 



(24) 



S+S*Sr = <4'(S+Sns;, (25) 
where form NN sequences, and 

S- 5* Sf = \ zz (Sr Sf) 5? • (26) 

We obtain 



M zz 

x\ z {u) = — - 2 — f—^ 

with M zz = ([i5 2 , Sl q ]) given by 

M zz = zC^- + (l - lq ) 

and 



(27) 



(28) 



«) 2 = J (1 " 7*){A" + 2zarc[° ) - + (l 7 ,)}, (29) 



A 22 = 2{5(5+l) - ((5 2 ) 2 ) 

+ [\ zz - {z + l)al z ]C { ^- + 

+ a?[(z-2)c[?- + + c!$-+]}. (30) 

As for the transverse correlations [cf. Eq. Ijlip], in the 
case S = 1/2 we have A 22 = 0. The correlation functions 
C^ zz are calculated fromi 7 - 



with 



M a 



■[l + 2n«)]. 



(32) 



Let us consider the magnetic susceptibility \ = ^^%Xs 
with xs = d(S z )/dh, which we denote by isothermal 
susceptibility, and its relation to the Kubo susceptibility 
(|27p. From the first and the second derivatives of the 
partition function with respect to h we obtain the exact 
relation 



Xs 



_ 1 M0)zz _ 1 ^ zz 

R 



(33) 



where C^' zz = C y ^' zz — (S z ) 2 , and the Fourier transform 



reads C zz 



(0)zz _ n (0)zz _ 

N(S z ) 2 5 q , . By Eqs. (27} and ((32} the 
uniform static Kubo susceptibility Xo 

may be expressed as Xq z 



ilimC 22 

1 9^0 9 



Urn lim Xq Z M 
= ilimC 22 = 

1 q^O q 



^Cq=o- That is, within our theory the isothermal and 
Kubo susceptibilities agree at arbitrary fields and tem- 
peratures. Using Eqs. l(2"7} to ((2"9"} we have 



d(S z 
dh 



2C 



(o)-+ 



(34) 



The equality (|34[) is an additional equation for determin- 
ing the parameters of the theory. 

Considering the ground state, at T = we have the 
exact results 



C 



(«)-+ 

R 



(0) 



0, c% )zz 



(0) 



5 2+,i . ((5 2 )™)(0) = S n . 

(35) 

The regularity conditions (jTTJ) read as go = h/zS. From 
go = 2 the field /io(0) is given by ho(0) = 2zS. Taking g 
from Eq. (fT5|) we get the equation A H = 2(1 — aj 1 )hS. 
This equation can be fulfilled only, if a{ (0) = 1 and 
A -1 (0) = 0, because in the ground state of the fer- 
romagnet at h ^ Q all quantities do not depend on h. 
Taking A H from Eq. (jTTJ) we get the parameter rela- 
tion A+-(0) + (z- l)c4~(0) = z - 1/2S. For 5 = 1/2 
(A -1 = 0) we have a\ (0) = 1. Concerning the zero- 
temperature values of a\ z and A 22 , they can be deter- 
mined only in the limit T — > 0, since Eqs. (|3ip and 
for C 



(0)Z2 

R 



contain M zz with limT^o M zz 



0. 



To evaluate the thermodynamic properties for arbi- 
trary spin, the transverse correlators C{q + > the longitu- 



z\n+l 



(n)zz \ 
10 

zz) have to be determined as 



C)'n >zz ), and the parameters 



dinal correlators (((S z ) 
and A^ (/ii/ 

solutions of a coupled system of self-consistency equa- 
tions for arbitrary temperatures and fields. Note that 
for S > 1/2 the parameters and A" M have not to be 
calculated separately, because they only appear in the 
combination given by A ufl . The correlation functions 
+ are calculated from the Green functions accord- 



q{o)zz _ J_ (jzz^iqR _|_ igz\2 ^2) ing to Eqs. (TT5} . To determine the 4(5+ 1) quantities 

V ^ q ^ ((5 2 ) n+1 ) and c[o )zz with n = oc 1 ~- J 



25 



1, aj , and 



5 



A" M , we have 65* + 3 equations, namely the regularity 
conditions (fTf]). the sum rules (HU) and (EU), Eqs. (JSJ) 

for ((S z ) 2 } and C$ )z *, and the equality (04]). That is, 
for S > 1/2 we have 25* — 1 more equations than quan- 
tities to be determined. To obtain a closed system of 
self-consistency equations for S > 1/2, i.e. to reduce the 
number of equations (in addition to those for C^, h ) 
to 4(S + 1), we consider two choices. First we take into 
account the higher-derivative sum rule (|21|) with n = 
only. As revealed by numerical evaluations, the specific 
heat of the ID model strongly deviates from the QMC 
data for S = 1, and for S > 1 it even becomes nega- 
tive at low fields and temperatures. Therefore, we adopt 
another choice, which yields a good agreement of all ther- 
modynamic quantities with the QMC data for S — 1 and 
which is used for S ^ 1 throughout the paper. Namely, 
we take into account the higher sum rules (|21|) with n = 

i(0)zz 



(i) S = ±: Using the identities (5f ) 2 = 1/4 and STf Sf = 
-iS7 (cf. Eq. ((12J>) the sum rules dUD and (J2T|) for 
n = simplify, where the higher sum rule (|2ip reduces 
to 



l(S z )~C w ] =-I^(-l)^+-A(%(o, gi ). (36) 



Note that this sum rule may be also obtained from the 
exact representation (|23p of the internal energy which in 
the case 5 = 1/2 becomes (cf. Ref. Gil) 

z h 

I /• + 00 7 

- nY,J ^(e q + u)Im((S+;SZ q )) u n(") (37) 



and with n = 1 instead of Eq. dSTJ for Cy" . To jus- with = ^ _ )/2 + ^ if ((5+;Sl )) w given by 
tify this choice within the theory itself, the correlator Eq ^ for n = Q ig inserted into Eq . The spec . 

C{ resulting from the closed system of equations is 



compared with c[^ zz calculated by Eq. (|3lj) . For ex- 
ample, in the ID S = 1 model at the fields h = 0.05 
and 0.1 the deviation is found to be less than 2% at all 
temperatures except for the region 0.1 < T < 1, where 
the maximal deviation is about 9% for T ~ 0.3 and 0.4, 
respectively. From the solution of the self-consistency 
equations in region I and from Eq. ([161) with go = 2 the 
boundary between regions I and II, ho(T), is determined. 
In Fig. [Q h (T) is plotted for S = 1/2 and S = 1. Note 
that in experiments realistic values of temperature and 
field lie in region I. Therefore, below nearly all results are 
presented in this region, and only some results for high 
enough temperatures and fields in region II are shown in 

Fig.m 

Let us finally make some comments on the evaluation 
of the theory for different spin values. 



tra u)+- and lj zz are given by Eqs. |10]), (TT]), ([29]), and 
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T 



(|30|) with A + ~ = and X zz = 0. We have to solve a 
closed system of coupled self-consistency equations for 



the seven quantities (S z ), C 



a/, and A v » (or 



i 2 L )- Note that in previous approache s i the simpli- 
fied choice ou^ — is taken disregarding the equality 
(|3~4"]) and not using either the condition (TT7]) (Ref. |T3) or 
the higher sum rule (J3SD (Ref. G3). 

(ii) 5^1: Let us specify the identities ([20]) and ([22]) 
which are used to reduce the sum rules (|19[) and ([2"T]) for 
n = 2S - 1, respectively, for S = 1 and S = 3/2. For 
S = 1 we have (S z ) 3 = Sf and (SffS^ = ~S Z S~, 
and for S = 3/2 we get (Sf) 4 = f(S'f) 2 - ^ and 

\3o- _ 



5(qz\2 9 
2V°i J 16 

(S z fSr = -l(S?) 2 S- + \S z Sr + §<?r. For 5 = i 
a closed system of coupled self-consistency equations for 



u 10 i a l ) 



and 



the ten quantities (S z ), ((S^) 2 ), C 10 
A 1 ^ has to be solved. 

In the case h = we have (5 Z ) = 0, and the correla- 
tors for n — only are needed. The spin-rotation sym- 
metry, implying + = 2C^ ZZ = Cr, is preserved 
by the second-order theory with — a i Z 2 = a i.2 an d 
A+- = X zz = A. Using ((S^) 2 ) =}S(S + 1) following 
from Eq. ©, the Eqs. ([10 ]) . ([IT ]) , (p f , and (|3CT ]l yield the 



spectrum lu~1 



uj„ given by 



^ = -(l-7 <? ){A + 2za 1 C 10 (l-7 q )} (38) 



with 



A = -S(S+l) + 2{\-(z + l) ai }C w 
+ 2a 2 {(z-2)C 11 + C 20 }, 



(39) 



FIG. 1: Boundary ho(T) in the h — T plane separating region 
I, h < ho(T), where the equality ui^ = h [cf. Eq. (|14[) ] may 
be fulfilled, from region II, ft > ho(T), where h > Ug~ for all 
<7- 



which agrees with the result of Ref.ll9l. if we put a 2 = oti. 

The susceptibility \s — Xo z resulting from Eq. ([34]) 
is given by \S = 2Cio/A. The correlators C R T)ZZ are 
calculated from Eqs. ([3"T]) and ([3^]) with (S z ) 2 replaced 
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by C zz = — Cr ZZ (see Refs. El, M), where the con- 
it 

densation part C zz describes long-range order (LRO). At 

T = we have the exact result C^f z = \S 2 . The ferro- 
magnetic LRO is reflected in the divergence of xs, so that 
A(0) = §5(5 + 1) + |5 ,2 {A -(z+ l)ai + (z- l)a 2 } = 
and uj q = z v /2ai/3S'(l - j q ). Then, by Eq. ([31]) we 

get Cft ZZ (0) = S/^/6ai5iifl + C zz resulting in the sum 
rule [R = 0, cf. Eq. ©] \s{S + 1) = S/y/6al + C zz , 
and in C 22 = ±5 2 (fl^O). Finally, at T = we obtain 
ai(0) = 3/2 and A(Q)+(z-l )a 2 (0) = ±(3z+l)-l/5. For 
5 = 1/2 (A = 0) we have a 2 (0) = ai(Q) = 3/2. At finite 
temperatures there is no LRO in the ID and 2D systems 
implying C zz = 0. The higher sum rule ([2"Tjl for n = or, 
equivalently, Eq. ([23]) turns out to be trivially fulfilled. 
Therefore, following Ref. we put a?2 = u\ = a and 
A(T) = A(0) = 2 - 1/5 and determine a(T) from the sum 

rule <^ 0)zz = ±5(5+1). 

III. QUANTUM MONTE CARLO 
SIMULATIONS 

In order to assess the accuracy of the approxima- 
tions employed in the Green-function theory presented 
in the previous section we perform QMC simulations. 
The Heisenberg ferromagnets with 5 = 1/2 and 5 = 1 
placed on chains or square lattices with periodic bound- 
ary conditions are simulated using the stochastic series 
expansion (SSE) method ) 27 ' 28 which utilizes the high- 
temperature series expansion 

00 on 

Z = Tr e-? H = EE^f( a l (-^)" I") ' ( 4 °) 

a n—0 

where the first sum is over a complete set of states | a) , 
usually taken as the eigenvectors of the 5 2 operator. 
By decomposing the Hamiltonian into diagonal and off- 
diagonal bond operators, introducing constant unit op- 
erators to assure positivity, and reexpanding (|40|) , one fi- 
nally ends up with a non-local loop representation which 
allows very efficient sampling] 27 ' 28 To minimize the ef- 
fect of self-crossing and back-tracking, the directed loop- 
updating scheme is employed. 

After initial thermalization with about 10 6 Monte 
Carlo steps, the measurements are made after each step. 
During the simulation, the energy, magnetization, and 
correlation functions are measured and stored in a time 
series file, from which the specific heat and magnetic 
susceptibility can be computed using the fluctuation- 
dissipation relation. The correlation lengths are ex- 
tracted from the exponential falloff of the correlation 
functions and for comparison also by means of the 
second-moment method. 29 Only for correlations smaller 
than one lattice spacing, small systematic deviations are 
visible. All those observables can be easily expressed by 
states of the spins on the lattice and the number and 



types of operators^ The whole simulation usually takes 
of the order of 10 7 Monte Carlo steps. The statistical 
error bars are estimated by the Jackknife method. 31 

The results presented in this paper are generated for 
5=1/2 chains of length up to L = 1024 and for 5 = 1 
up to L — 64. In two dimensions we simulate square 
lattices of edge length up to L — 64. By comparing 
the results for different lattice sizes we made sure that 
for the investigated range of temperatures and fields, the 
thermodynamic limit of the considered observables lies 
within the statistical error bars of the numerical results. 



IV. RESULTS 

As described in Sec. [Til the quantities of the Green- 
function theory determining the thermodynamic prop- 
erties have to be calculated numerically as solutions of 
a coupled system of non-linear algebraic self-consistency 
equations. To this end, we use Broyden's method 3 -? which 
yields the solutions with a relative error of about 10~ 7 
on the average, where the numerical error increases with 
decreasing field and temperature. The momentum inte- 
grals occurring in the self-consistency equations are done 
by Gaussian integration. Considering the 5=1/2 fer- 
romagnet, in Refs. [l?] and [3 the thermodynamic quan- 
tities, except for the transverse and longitudinal correla- 
tion lengths, are calculated. Therefore, we present only 
some results for 5 = 1/2 (see Figs. [31 and [T5")) which 
visibly improve those of Ref. Il7l 



A. Magnetic susceptibility 

Let us first consider the susceptibility xs m the case 

h = 0, xs — 2Cio/A ( see Sec. [ITJ). In one dimension, 

2 

the low-temperature expansion yields lim xsT 2 = — 5 4 

(Ref. GJ). Note that this result agrees with that ob- 
tained by the modified spin- wave theory (MSWT). 33 
For 5 = 1/2 we have hrnxsT 2 = 0.041667 which is 

in very good agreement with the Bethe-ansatz value 
hinxsT 2 = 0.041675 (Ref. Hi). On the other hand, pre- 
vious QMC simulations by Handscomb's method on an 
N = 256 chain combined with a renormalization-group 
approach 3 ^ yield hm xsT 2 = 0.0329 (note that x plot- 
ted in Ref. [H and defined in Ref. [36| is related to xs by 
X = 3xs/5 2 ). To resolve the discrepancy between the 
QMC results of Ref. |35j and the Bethe-ansatz value, we 
perform QMC simulations for chains up to N = 1024 
sites. The results at very low temperatures are shown 
in Fig. [2] (taking the same plot as in Ref. l35h and com- 
pared with the Bethe-ansatz data^ the QMC data of 
Ref. HU, and with the Green-function theory. Above a 
characteristic temperature, which decreases with increas- 
ing chain length, our QMC data agree very well with the 
Bethe-ansatz results. On the contrary, the QMC results 
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FIG. 2: Zero-field susceptibility of the ID 5=1/2 ferromag- 
net. The results of the Green-function theory (solid line) and 
the QMC data (+ : L = 256, x : L = 1024) are compared 
with the QMC data of Ref. [35j (■) and the Bethe-ansatz re- 
sults of Ref. [34| (°). In the inset the finite-size scaling of 
the zero-temperature limit of xsT 2 calculated by QMC is de- 
picted. The dashed line shows the least-square fit of the data 
by a linear dependence. 



of Ref. 35 for XsT 2 are lower than ours by 4% on the av- 
erage. To determine the limit lim xsT 2 from our QMC 

data, we perform a finite-size scaling analysis. To this 
end, for each chain length we linearly extrapolate the 
low-temperature linear part of the curve xsT 2 to T = 
and fit the limiting values as function of 1 /N by a linear 
dependence (see inset of Fig. [2]). The extrapolation to 




0.1 0.2 0.3 

T 

FIG. 3: Susceptibility of the ID 5 = 1/2 ferromagnet at 
h — 0.005 and 0.05, from top to bottom, where the results 
of the Green-function theory (solid lines) and of the Green- 
function method of Ref. [l8| (dashed lines), the QMC data 
(filled symbols, L = 128), and the Bethe-ansatz results of 
Ref. [13 (open symbols) are shown. In the inset the ID mag- 
netization at ft = 0.005 and 0.05, from bottom to top, is 
depicted. 
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FIG. 4: Magnetization of the (a) ID and (b) 2D S = 
1 ferromagnet in magnetic fields of strengths (a) h — 
0.1,0.2,0.4,0.6,1.0, and 2.0, from bottom to top and (b) 
h = 0.005,0.01,0.05,0.1,0.5, and 1.0, from bottom to top, 
as obtained by the Green-function theory (solid lines) and 
the QMC method for L = 64 (•), compared with RPA results 
(dashed lines). 



l/N = yields lim lim y s T 2 = 0.0413 ± 0.0005 which 

agrees, within the given statistical error, with the Bethe- 
ansatz value. 

The 2D zero-field susceptibility in the second-order 
Green-function theory increases exponentially for T — > 
0, xs oc exp(2itS 2 /T) (Ref. 0), where the exponent 
is smaller by a factor of two as compared with that 
found in the MSWT22, and in the renormalization-group 
approach^ 

Now we consider nonzero fields and calculate the sus- 
ceptibility xs — d(S z ) jdh. First we show the magnetiza- 
tion. For 5 = 1/2, as an example, (S z ) in the ID model 
is depicted in the inset of Fig. [3] For the S — 1 ferromag- 
net our analytical and QMC results in comparison with 
the RPA are plotted in Fig. [5] Let us emphasize the ex- 
cellent agreement of the theory for the chain (Fig. SJa)) 
with the QMC data over the whole temperature and field 
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FIG. 5: Susceptibility of the ID S — 1 ferromagnet (a) at low 
fields, h = 0.005, 0.01, 0.03, and 0.05, from top to bottom, and 
(b) at higher fields, h = 0.1, 0.2, 0.4, 0.6, 1.0, and 2.0, from top 
to bottom, where the Green- function (solid lines), the QMC 
(•, L = 64), and the RPA results (dashed lines) are shown. 



FIG. 6: Susceptibility of the 2D S — 1 ferromagnet (a) at 
very low fields, h — 0.005 and 0.01, from top to bottom, and 
(b) at higher fields, h — 0.05,0.1,0.5, and 1.0, from top to 
bottom, obtained by the Green-function theory (solid lines), 
QMC for L = 64 (filled symbols), and RPA (dashed lines). 



regions. For the ID ferromagnet the RPA is a remark- 
ably good approximation for (S 12 ), as was also found in 
the case S = l/2»i£ In two dimensions (Fig. SJb)), as 
compared with the QMC data, the results of our theory 
at higher temperatures are somewhat worse than those 
of the RPA. This is in contrast to the 2D S = 1/2 ferro- 
magnet for which we obtain slightly better results than 
the RPA at all temperatures and fields (improving our 
previous findings^). 

The susceptibility for h ^ vanishes at T — 0. There- 
fore, Xs(T) has a maximum at where increases 
and the height of the susceptibility maximum XsiT^) de- 
creases with increasing field. For S = 1/2, in Fig. [3] the 
low-field susceptibility in the ID model is shown, where 
for h = 0.005 a better agreement of the theory with the 
Bethe-ansatz results is found than in Ref. [ItI Note that 
our QMC data are in a very good agreement with the 
Bethe results. For comparison, in Fig. [3] the susceptibil- 
ity in the simplified approach with ot^ = a\ M (Ref. [H), 
where the equality (I34[) is disregarded and the regular- 



ity condition (|17p is used instead of the higher sum rule 
([55]) . is plotted as well. It is remarkable that xs in this 
approach is in a better agreement with the exact meth- 
ods than the susceptibility in our extended theory with 
a^ 1 ^ a^. However, considering the correlation length 
the situation changes qualitatively (see below). For S = 1 
the susceptibility is plotted in Figs. and [5] In one di- 
mension (Fig. 0, the good agreement between Green- 
function theory and QMC corresponds to the results de- 
picted in Fig.QJa). As compared with the QMC data for 
the 2D model (Fig. in RPA the maximum position 
is somewhat better reproduced than in our theory. 
To analyze the field dependence of TJ* and xs(T^) in 
more detail as in our previous paper^i the calculations 
are extended to a much broader field region, 0.001 ^ h ^ 
10. As can be seen in Fig. [T^a), at low fields the theory 
may be well fitted by the power law 

T* = ah\ (41) 
where the field regions and the values of a and 7 are 
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TABLE I: Validity regions (h) and coefficients of the power 
laws (|4"Tj) and (TJ2| for the susceptibility of the ID and 2D 
S = 1/2 and S = 1 ferromagnets. 





S = 


1/2 


S = 


1 




ID 


2D 


ID 


2D 




0.001 - 1.0 


0.001 - 0.1 


0.001 - 2.0 


0.001-0.1 


a 


1.013 


1.149 


1.823 


2.433 


7 


0.596 


0.192 


0.565 


0.144 


h 


1.0-10.0 


1.0 - 10.0 


2.0-10.0 


2.0-10.0 


a 


0.661 


0.666 


0.917 


0.929 


b 


0.443 


0.961 


1.136 


2.494 




0.001 ~~ 0.01 0.1 1 10 

h 

FIG. 7: Field dependence of the (a) position and (b) height of 
the susceptibility maximum obtained by the Green-function 
theory for the 5 = 1/2 (•) and S = 1 (o) ferromagnets and fit 
by power laws (solid lines) in comparison with the QMC data 
(+, S = 1, L = 64). The inset shows the fit of T* at high 
fields by a linear dependence. For clarity, xs(T^) is plotted 
for S = 1 only. 

given in Table HI Let us point out that the theory for 
the ID S — 1/2 model is in reasonable agreement with 
the Bethe-ansatz result at h ^ 0.1, 17 a — 0.765 and 
7 = 0.576. In the high-field region, obeys a linear 
dependence (cf. inset of Fig. Ufa)), 

TX=ah + b (42) 

with a and b given in Table |TJ Note that the linear law 
(|42|) was not found in Ref. [l?]- Our results for the max- 
imum height Xs(^m) as a function of h may be well de- 
scribed in the whole field region 0.001 ^ h ^ 10.0 (see 
Fig. Efb)) by the power law 

X s(TX) = bh ^ (43) 

where the coefficients are given in Table [TTJ The values 
of b and j3 for S = 1/2 slightly deviate (by about 5% 
on the average) from those found previously^ Again, 



TABLE II: Coefficients of the power law (143 [) for the suscep- 
tibility of the ID and 2D S = 1/2 and S = 1 ferromagnets in 
the field region 0.001 ^ h «S 10.0 





S = 


1/2 


S = 


1 




ID 


2D 


ID 


2D 


b 


0.192 


0.166 


0.362 


0.305 


P 


-0.925 


-0.850 


-0.941 


-0.867 



our theory for S = 1/2 is in reasonable agreement with 
the ID Bethe-ansatz result at h ^ 0.1, b = 0.208 and 
13 = -0.952 (Ref. 03). 

For comparison, we consider the power-law behavior 
in RPA. We find the RPA results in the low- and high- 
field regions to be well fitted by the laws PT^ - PS")) . where 
the coefficients are in good agreement with the values 
given in Tables [T] and [TTJ More precisely, for the ID 
and 2D S = 1/2 and S = 1 models the average de- 
viations of the coefficients in the laws ([4"Tj) . (|4"2"]) , and 
(U2J) amount to about 6%, 3%, and 2%, respectively. For 
example, considering the 5=1/2 ferromagnet in high 
fields, 2 ^ h ^ 10, we obtain the linear dependence (l4"2"|) 
for the ID (2D) case with 5, = 0.657 (0.661) and b = 0.496 
(1.015) which yields a better fit than the power law (|41[) . 
Recently, in Ref. l37| such a law was given for the ID (2D) 
model in the region 3 (4.4) ^ h ^ 6.5. Even in this lim- 
ited field region, we find the fit by the linear law (|42]) to 
be slightly better than the fit by the power law (|41 [) (sec 
Ref. [33). 



B. Correlation length 

To obtain the transverse and longitudinal correlation 
lengths ^ and £ zz we consider the long-distance cor- 
relators C^ 0) ~ + and C [ ° )zz = C [ £ )zz - (S z ) 2 with C ( ° )zz 
calculated by Eq. (|3ip . respectively. Note that the tem- 
perature dependence of both C^' + and C < j^ >zz exhibits 
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a maximum, because the correlators vanish at T = 0, 
following from Eqs. (f3"Tj) and (f3"5)) . and for T — > oo. By 
the asymptotic ansatz 



C 



(o)-+ _ 4 +- 



R 



= A+- exp(-R/?~), 



Ft. 



A-exp(- J R/r z ), 



(44) 



(45) 



and the logarithmic plot of the correlators as functions 
of R = \R\ the inverse correlation lengths are evaluated 
numerically from linear fits. 

In the literature, often the correlation length is deter- 
mined from the expansion of the static spin susceptibility 
around the magnetic wavevector (see, e.g., Refs. [HI 
andl20l). In the ferromagnetic case we expand the static 
susceptibilities Xq~ (resulting from Eqs. JTOJ) - (JT3J) , (fT4")) . 



and (jT5j) ) and \ 



Aq 



and 



xo7[i 



(given by Eqs. 
(vfj, = 



(£^)V 



zz). 



<Ja+-(S*)/h 



2af c[° > ~ + /A z 



around q = 0, 
We obtain 

(46) 



(47) 



Deriving Eq. (|4"6"|) the regularity condition (jTTJ) for n = 0, 
which reads as /i(S' z ) = zCiogo, and Eq. (fl"6)l . yielding 
the relation A+~ = 2/i(Ci /(S z ) - a+~(S z }), have been 
used. Let us point out that the correlation lengths £J^ 
generally deviate from £^ M defined by Eqs. (I44[) and (|4"5f . 

First we consider the correlation length in zero field, 
where £ H = £ zz = £. In one dimension, the low- 
temperature expansion yields liniT->o £,T = S 2 (Ref. 
which agrees with the MSWT result 3 - 3 - and, for S — 
1/2, with the result obtained by the thermal Bethe- 
ansatz method of Ref. l38l . The renormalization-group ap- 
proach of Ref. [H combined with QMC simulations yields 
\im T ^ Q £T = 1.14S* 2 . In Fig. [5] the zero-field correlation 
length of the ID fcrromagnet is shown. Let us stress the 
very good agreement of our QMC data for S = 1/2 with 
the Bethe-ansatz results of Ref. |38|. Even on the finer 
scale of the inset, deviations are almost invisible. For 
comparison, also the QMC data of Ref. 35 and a one- 
parameter fit are given in the inset. Moreover, we obtain 
a good agreement of the Green-function theory, where 
£ is calculated from the definition (j4~4"|) , with our QMC 
data. In addition to £, in Fig. [8] the correlation length £ x 
calculated for S = 1/2 and S = 1 by Eq. gT]) [a zz = a, 



C 



(o)-+ 



in 



Cio, A zz = A given by Eq. ((MD] is plotted. 
For T < 0.25, i.e. £ > 1, £ x nearly coincides with £. 
With increasing temperature, i.e., with decreasing £ < 1, 
the deviation of £ x from £ appreciably increases. In the 
high-temperature limit we get £^ 1 = {3T/S(S + 1)} 1/2 
resulting from C w = 2{S(S + 1)] 2 /9T (Ref. IH). In the 
following we plot £ x in such cases only, where £ x remark- 
ably deviates from £. 
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FIG. 8: Zero-field correlation length of the ID ferromagnet 
with 5 = 1/2 obtained by the Green- function theory (solid 
lines) and by QMC simulations (x, L = 32) and with S — 1 
resulting from the theory (long-dashed line). For comparison, 
the correlation length £ x determined from the expansion of 
the static susceptibility around q — is plotted for S — 1/2 
(dotdashed line) and S = 1 (dotted line). The results for 
S = 1/2 are compared with the Bethe-ansatz data (o) of 
Ref. |M and with the QMC data of Ref.Oj! (•) depicted in the 
inset together with a one-parameter fit (short-dashed line). 



In two dimensions, the zero-field correlation length in 
the second-order Green-function theory increases expo- 
nentially for T -> 0, £ cx exp(7r,S' 2 /T) (Ref. [l]J. As is 
the case for the magnetic susceptibility, the exponent is 
smaller by a factor of two as compared with the MSWT 33 
and the renormalization-group approach. 36 

For h ^ the transverse and longitudinal correlation 
lengths reveal qualitatively different temperature depen- 
dences. Considering the transverse correlation length 
£^ shown in Fig. [9) the magnetic field cuts off the diver- 
gence of the zero-field correlation length at T — which 
corresponds to the absence of a phase transition and is 
evident from Eq. (|46j), £+ ~(T = 0) = ^/Sjh agreeing 
with the RPA result (|A.4j) derived in the Appendix. As 
can be seen in the inset of Fig. [HJa), in the ID 5=1/2 
model we obtain a good agreement of our analytical re- 
sults for T = 0.4 and h ^ 1.2 with the Bethe-ansatz data 
of Ref. [H. However, the comparison of the theory with 
the available Bethe data for T = 0.4 and fields up to 
h = 4 and for T = 0.2 (Ref. l39t ) is hampered by numeri- 
cal uncertainties resulting from too small values of A^ . 
Note the remarkably good agreement of £^ with the 
RPA results (see inset). Concerning the dimensional de- 
pendence, in contrast to the case h — 0, £^ in one and 
two dimensions exhibits qualitatively the same behavior 
as T -> 0. In the 2D model [Fig. Mh)}, the deviation 
of £ x ~ from £^ increases with decreasing temperature, 
i.e., with increasing £^ > 1 which is clearly seen at 
h = 0.01 and is opposite to the behavior in the h = 
case. 

In Fig. [TU] the longitudinal correlation length of the ID 
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FIG. 9: Transverse correlation length of the (a) ID and (b) 
2D ferromagnet with S =1/2 (solid lines) and S = 1 (dashed 
lines) in the fields h = 0.01 and 0.1, from top to bottom. In 
the 2D case at h = 0.01, the correlation length calculated 
from the static susceptibility is shown for S — 1/2 (dotdashed 
line) and S = 1 (dotted line). In the inset the results of the 
Green-function theory are compared with the Bethe-ansatz 
data of Ref. H<3 (o) and the RPA (dotted line). 



FIG. 10: Longitudinal correlation length of the ID ferromag- 
net with (a) S = 1/2 and (b) S = 1 in the fields (a) h = 0.005 
and 0.05 and (b) h = 0.05 and 0.1, from top to bottom, calcu- 
lated by the Green-function (solid lines) and QMC methods 
(x, o; L = 32 and ; +; L = 32) and, for S = 1/2, by the 
method of Ref. Qjl (dotdashed lines). The inset exhibits the 
results for S = 1/2 at the strong fields h = 1, 3, and 5, from 
top to bottom, in comparison with the correlation length 
(dashed lines) obtained from the static susceptibility. 



ferromagnet is shown, where the QMC data are found 
to be in a fair agreement with our theory. This refers, 
in particular, to the 5=1/2 model, where our results 
obtained by the simplified approach of Ref. \±B are plot- 
ted as well. Considering h = 0.05, at low tempera- 
tures those results remarkably deviate from the QMC 
data and our extended theory with ^ a-^ . In con- 
trast to £ H , the behavior of £ zz as T — > is not con- 
clusive which is due to numerical uncertainties at low 
temperatures, where the long-distance correlators (To 
needed to calculate £ zz are very small. For example, 
for S — 1/2 and strong fields [see inset of Fig. ITDTa)] 
the relevant correlators in the temperature region, where 
results are not given, are smaller than about 10 -10 to 
10 -14 . Moreover, for S = 1 the results of the theory 
are reliable only at T > T ~ 0.1 and 0.3 for h = 0.05 
and 0.1, respectively [see Fig. [TDTb)]. At T < T , the 



relevant correlators, being smaller than about 10~ 4 , re- 
veal an unreasonable behavior. This may be ascribed to 
our choice of a closed system of self-consistency equa- 
tions for S > 1/2, as described in Sec. [Til Whereas the 
relative deviation of the NN correlators C 10 resulting 
from the self-consistency equations and from Eq. (f3"Tj) is 
small (see Sec. HTJ, the corresponding deviation of the 
correlators C^ zz becomes very large at low tempera- 
tures. Depending on the field and spin, the temperature 
dependence of £ zz in the ID ferromagnet reveals a max- 
imum at Tjj > 0. This anomaly can be clearly seen in 
the ID S = 1 model at low fields [Fig. [TTj^b)]. On the 
other hand, in the ID S =1/2 model the maximum ap- 
pears at high fields, h > 0.8 [see inset of Fig. [TUTa)]. 
Moreover, as can be seen from Fig. [TUl keeping the field 
h = 0.05 fixed, the maximum develops with increasing 
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spin. Note that a maximum of £ zz at a finite tempera- 
ture is not obtained by the approach of Ref. To our 
knowledge, such an anomaly in the correlation length 
has not been found before. To get some insight into the 
maximum of £ zz , we first suggest that larger correlation 
lengths may be connected with larger correlation func- 
tions. Correspondingly, we consider the maximum of 
C^ )zz at T^{R), where T^{R) > T^. By a detailed 
analysis we find T^(R) in the limit R — ► oo to coin- 
cide with Tjj in all cases, where £ zz has a maximum 
at T« > (see Fig. [TOD, i.e., lim*^ T zz (i?) = l£. 
This result is corroborated by the conditions for a max- 
imum which may be derived from the ansatz (|45|) . We 
get ±dhxC$ zz /dT = ±d\nA zz /dT + ±01nf/dT. At 
T™(R) we have ±dln£/dT = -±dk\A z * /&T and, for 
R — > oo, d^/dT — 0. As can be easily verified, the maxi- 
mum condition d 2 C ( n )zz /dT 2 < results in d 2 £ zz /dT 2 < 
0. To compare the QMC and Green- function methods 
yielding the anomaly of £ zz in the ID S = 1 model 
[Fig.fTOfb)] in more detail, in Fig.[TT]the distance depen- 
dence of the corresponding correlator C^ zz at h = 0.05 
is depicted. For T — 0.5 a very good agreement of both 
methods is found. 

In two dimensions, the anomaly of £ zz in the 5=1/2 
ferromagnet is more pronounced than in the ID sys- 
tem and appears already at low fields, as can be seen 
in Fig. [TO] In contrast to the ID case, both the QMC 
data and the Green-function theory clearly reveal a min- 
imum in addition to the maximum. Note that the sta- 
tistical QMC errors in the interesting temperature re- 
gion are smaller than the size of the symbols. Figure [TO] 
demonstrates the qualitative effects of our extended the- 
ory (a^ oty 1 ) on the temperature dependence of £ zz 
as compared with the simplified approach (a^ 
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FIG. 11: Correlation function C^ )zz = (5g5|j) - <5 Z ) 2 vs 
R = \R\ for the ID 5 = 1 ferromagnet in the field h = 0.05 
at T = 0.5 and 1.0, from top to bottom, calculated by the 
Green-function theory (open symbols) and by QMC (filled 
symbols, L = 32). 




FIG. 12: Longitudinal correlation length of the 2D 5 = 1/2 
ferromagnet at h — 0.05 calculated by the Green-function 
theory (solid lines), QMC simulations (•, L — 16), and by 
the method of Ref. LL8| (dashed lines) . In the inset the corre- 
sponding magnetization is plotted. 



Whereas this approach yields a slightly better agreement 
of the magnetization with the QMC data (see inset), it 
fails to describe the minimum-maximum anomaly. 

Figure [TO] shows the field and spin dependence of the 
temperature behavior of £ zz in the 2D ferromagnet. As 
results from the theory, the anomaly of £ zz becomes more 
pronounced with decreasing field and with increasing 
spin. Let us point out that our QMC data for h = 0.05 
yield a minimum and a maximum of £ zz for both the 
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T/S(S+1) 
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FIG. 13: Longitudinal correlation length of the 2D ferromag- 
net with 5 = 1/2 (solid lines) and 5=1 (dashed lines) in 
the fields h = 0.01, 0.05, and 0.10, from top to bottom, as 
compared with the QMC results at h = 0.05 for 5 = 1/2 (•, 
L = 16) and 5=1 (+, L = 16). The inset shows the correla- 
tion function C ( n )zz = (5o5«) - (S z ) 2 vs R = \R\ for the 2D 
5=1/2 ferromagnet in the field h = 0.05 at T = 0.4 and 0.6, 
from bottom to top, calculated by the Green-function theory 
(open symbols) and by QMC (filled symbols, L = 16). 
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S = 1/2 and S = 1 models and give confidence in the 
results of the theory. As in the ID model, the maxi- 



(0)2 

R 



mum of (f z at is related to the maximum of C 
by limft^oo T^ z (i?) = in all cases shown in Fig. Q]|] 
The minimum of £ 2Z results from the different tempera- 
ture dependences of C^) zz and A zz in the ansatz (|4"S")) . 
In analogy to Fig. [TTJ for a more detailed comparison, 
the inset exhibits the correlator C^ zz for 5 = 1/2 and 
h = 0.05 as function of the distance. The relative magni- 
tude of the correlators at T — 0.4 and 0.6 may be under- 
stood by the maximum in the temperature dependence 



of C 



(0)zz 
R 



C. Specific heat 



Let us first consider the NN spin correlation func- 
tions C^o + and C^ zz entering the internal energy 
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h(S z ). As an example, for 



the ID S = 1 model they are depicted in Fig. rjj] where 
we obtain a very good agreement of the analytical results 
with the QMC data. On the contrary, the RPA results for 
C 



(0)- 
10 



+ remarkably exceed the QMC data, and for c[^ zz 
the RPA yields negative values being incompatible with 
the ferromagnetic SRO. 

In Fig.[l5]the specific heat C = du/dT for the ID S = 
1/2 ferromagnet at low fields is plotted. Again, our QMC 
data agree very well with the Bethe-ansatz results^ At 
very low magnetic fields, the low-temperature maximum 
appearing, in the exact approaches at h < 0.008, in ad- 
dition to the high-temperature maximum is much better 
described by the theory than we have found in Ref. fl7l . 
In our Green-function theory this maximum appears up 
to higher fields, h ^ 0.071, and the deviation of the max- 
imum position l from the Bethe-ansatz and QMC 
values in the region 0.001 ^ h ^ 0.01 is less than 8%. 
Considering very low fields, h = 0.001 to 0.01 in steps of 
0.001, Tf n l and the height C(T% A ) are fit by the power 
laws 



0.462 ti 3 



CCZ^i) = 0.394 h° 



(48) 



The exponents are in good agreement with the values 
of the Bethe-ansatz results,! 7 . T% A = 0.596 ft, 542 and 
C(T° A ) = 0.513 h°- 22S . Note that the specific heat in 
the 2D model has only one maximum^ 

Figure [TBI displays the specific heat of the ID 5 = 1 fer- 
romagnet. At low magnetic fields, 0.007 < h < 0.057, be- 
sides the high-temperature maximum, a low-temperature 
maximum appears (see Fig. |16f a)). The position 1 
of this maximum obtained by the Green-function theory 
nearly agrees with the QMC results. As in the 5=1/2 
case) 1 - 7 in RPA a double maximum is not obtained (see 
inset of Fig. rTBT a)). and the values of the specific heat 





FIG. 14: (a) Transverse and (b) longitudinal nearest-neighbor 
two-spin correlation functions of the ID S = 1 ferromagnet at 
the fields h = 0.1, 0.2, 0.4, 0.6, 1.0, and 2.0, from left to right, 
obtained by the Green-function theory (solid lines), QMC (•, 
L = 64), and RPA (dashed lines). 



FIG. 15: Specific heat of the ID S = 1/2 ferromagnet ob- 
tained by the Green-function (solid lines) and QMC (filled 
symbols, L = 128) methods at low fields, h = 0.005, 0.03, and 
0.1, from bottom to top, and compared with the Bethe-ansatz 
data of Ref. Il7| (open symbols). 
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FIG. 16: Specific heat of the ID S = 1 ferromagnet obtained 
by the Green-function theory (solid lines) (a) at low fields, 
h = 0.005,0.01,0.03, and 0.05, from bottom to top, with the 
QMC results for L = 64 (filled symbols) and (b) at higher 
fields, h = 0.1,0.2,0.4,0.6, 1.0, and 2.0, from left to right, in 
comparison with the QMC results for L — 64 (filled symbols) . 
The inset shows the RPA data at the fields given in (a), from 
top to bottom at T = 0.1. 



maximum are much higher than the QMC values which 
is ascribed to a poor description of SRO in RPA (see also 
Fig. [H|). The specific heat of the ID 5 = 3/2 ferromagnet 
is shown in Fig. 1171 There is no low-temperature maxi- 
mum, but only a hump at low enough fields. For higher 
spins qualitatively the same behavior is found. The spe- 
cific heat for the 2D 5 = 1 ferromagnet is plotted in 
Fig. [H] As in the case 5 = 1/2 p in two dimensions 
only one maximum appears. At small fields the position 
of the maximum in the Green-function theory is remark- 
ably shifted to higher temperatures as compared with 
the QMC data. Note that the RPA curves at low fields 
(see upper inset of Fig. [T8]) exhibit a too large maxi- 
mum height, as was also found in the ID model (inset of 
Fig. (Ma)). 

From our investigations of the maximum behavior of 
the specific heat in dependence on spin and dimension we 



FIG. 17: Specific heat of the ID S = 3/2 ferromagnet calcu- 
lated by the Green-function theory at h = 0.01, 0.1 and 1.0, 
at T = 1.5 from bottom to top. 



conclude that the appearance of two maxima is a distinc- 
tive effect of quantum fluctuations which decrease with 
increasing spin and dimension. Note that in ferromag- 
nets quantum fluctuations occur at nonzero temperatures 
only, whereas in antiferromagnets they are important al- 
ready at T = 0. The characterization of the occurrence 
of two maxima in the temperature dependence of the 
specific heat of the Heisenberg ferromagnet as a peculiar 
quantum effect is corroborated by recent QMC simula- 
tions of the ID classical Heisenberg model and the ID 
5=1/2 Ising model in a magnetic field, 40 where only 
one maximum in the specific heat was found. 




FIG. 18: Specific heat of the 2D 5=1 ferromagnet at h = 
0.01 and 0.05, from bottom to top, and, as depicted in the 
lower inset, at h = 0.1 and 1.0, from left to right, where the 
Green-function (solid lines) and QMC (filled symbols, L = 64) 
results are shown. In the upper inset the RPA results for 
h — 0.01 and 0.05, from top to bottom, are plotted. 
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D. Comparison with experiments 

Let us compare our results with experiments on 5 = 
1/2 quasi-lD ferromagnets, where we focus on the pos- 
sible observation of two maxima in the temperature de- 
pendence of the specific heat as a characteristic feature 
of ID ferromagnets in a magnetic field. 

The copper salt TMCuC [(CH 3 ) 4 NCuCl 3 ] was 
shown 2 .^ to be a good ID Heisenberg ferromagnet which 
is reflected in the small value of the Neel temperature 
Tyv = 1.24K for 3D ordering. 3 Determining the exchange 
energy J by a least-squares fit of the theory for 5 = 1/2 
to the experimental data for the magnetization as a func- 
tion of the magnetic field H at T = 4.1K^ we obtain 
J = 6.18meV and a very good agreement with exper- 
iments, as can be seen in the inset of Fig. [T9j Note 
that the value of J lies between the values given in 
Ref. (J = 5.17meV) and in Ref. 1 (J = 7.76meV). 
According to the QMC and Bethe-ansatz results for 
the ID 5 = 1/2 ferromagnet, two maxima of the spe- 
cific heat occur for h < 0.008 or, using the relation 
h = 1.16 x 10- 2 iJ[kOe]/J[meV], for H ^ 4kOe. In 
Fig. [19] the specific heat, as predicted by the theory us- 
ing the fit value of J, is plotted. The low-temperature 
maximum for H = 2kOe, 3kOe, and 4kOe occurs at 
T% A = 2.0K, 2.5K, and 2.9K, respectively. The high- 
temperature maximum (not shown in Fig. I19p appears 
at about Tg 2 = 37.4K with C{Tg 2 ) = 1.18J/molK for 
all fields considered. In the quasi-lD system the anomaly 
of the specific heat at TV, which cannot be described by 
our theory for a purely ID system, may mask the low- 
temperature maximum, if 1 is not sufficiently larger 
than T N . At H = 3kOe (4kOe) we have T^ JT N = 2.0 




FIG. 19: Specific heat of the copper salt TMCuC (Refs. @ 
and 0, Neel temperature Tm — 1.24K), as predicted by the 
theory for the S — 1/2 ID ferromagnet in the magnetic fields 
H = 2kOe, 3kOe, and 4kOe, from bottom to top, with J = 
6.18 meV obtained from the fit of the reduced magnetization 
m = m(H)/m{H = 8.7kOe) at T = 4. IK to experimental 
data (o) shown in the inset. 



(2.3). From this we predict that in TMCuC above T N 
two maxima in the specific heat at moderate magnetic 
fields, H = 3 — 4kOe, may be observed. 

Considering the quasi-lD organic ferromagnet p- 
NPNN (Ci 3 Hi 6 N 3 04) in the 7 phase with J = 
0.37meV^£ where the phase transition at T/v — 0.65K 
for H = persists up to H = 1.8kOe (Tjy ~ 0.5K), 
two maxima of the specific heat above T/v cannot be ob- 
served, because, at ft, < 0.008, (H < 0.26kOe), we have 
1 ~ 0.19K < Tn. The analogous situation, in which 
the low-temperature maximum in the specific heat of the 
ID ferromagnet cannot be seen, is found for the following 
compounds. Considering the ferromagnetic chains in the 
quasi-lD magnet /5-BBDTA GaBr 4 with J = 0.375meV£ 
we have T^ ll < 0.19K which is lower than the temper- 
ature of the specific-heat cusp, Tq 5? 0.4K, caused by 
the interchain coupling. For the CuCbj-TMSO (tetram- 
ethylsulfoxide) [DMSO (dimethylsulfoxide)] salts with 
J = 3.36 [3.88]meVl we get T% tl < 1.7 [1.96]K being 
lower than the temperature of the susceptibility maxi- 
mum, 3.9 [5.4]K, indicating the influence of the antifcr- 
romagnetic interchain coupling. 



V. SUMMARY 

In this paper we have developed a second-order Green- 
function theory for the ID and 2D Heisenberg ferro- 
magnets in a magnetic field which extends our previous 
approach^! to arbitrary spins and by the calculation of 
the correlation length. In addition, we have performed 
QMC simulations of the 5 = 1/2 and 5=1 models on a 
chain up to N — 1024 sites and on a square lattice up to 
TV = 64 x 64 using the stochastic series expansion method 
with directed loop updates. The approximate analyti- 
cal and quasi-exact numerical results turned out to be 
in good agreement, in particular for the ferromagnetic 
quantum spin chains. Analyzing the field dependence 
of the maximum in the temperature dependence of the 
magnetic susceptibility over a much broader field region 
as considered previously^ we have found power laws for 
the position and height of the susceptibility maximum. 
The transverse and longitudinal correlation lengths were 
shown to have qualitatively different temperature depen- 
dences. Depending on spin, field, and dimension, the 
longitudinal correlation length ^ zz reveals an unexpected 
anomaly: with increasing temperature, £ zz exhibits a 
minimum followed by a maximum. By a detailed inves- 
tigation of the specific heat of the Heisenberg chain with 
arbitrary spin, two maxima in its temperature depen- 
dence at low magnetic fields were detected for 5=1/2 
and 5 = 1, whereas for 5 > 1 only one maximum ap- 
pears, as in the 2D case. The existence of two specific- 
heat maxima was identified as a distinctive quantum ef- 
fect. The theory was compared with magnetization ex- 
periments on the ID copper salt TMCuC, and predictions 
for the temperature dependence of the specific heat, in 
particular for the occurrence of two maxima, were made 
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which should be measurable experimentally. 

AKNOWLEDGMENTS 

The authors wish to thank J. Richter, N. M. Plakida 
and S. Wenzel for valuable discussions. This work was 
partially supported (L. B. and W. J.) by the EU through 
the Marie Curie Host Development Fellowship under 
Grant No. IHP-HPMD-CT-2001-00108 and the super- 
computer time grant No. hlzl2 of the John von Neu- 
mann Institute for Computing (NIC), Forschungszen- 
trum Julich. 



APPENDIX: RANDOM-PHASE 
APPROXIMATION 

It is of interest to compare our results for finite mag- 
netic fields with the RPAi^S. Considering the equation of 

motion |(3J) the Tyablikov decoupling iS q = o-> q Sq yields 
SM-)) U = , u; q = z(S*)(l - lq ) + h, 

y y bJ — UJ q 

(A.l) 

with M( n ) + given by Eq. ((4]). Comparing the correla- 
tion function ((S z ) n S~ S+) resulting from Eq. (jA.ip with 
the expression obtained by Eq. © multiplied by (Sf ) n 

s 

and using the identity JJ_ (Sf - m) = 0, (S z ) is ob- 

m=-S 

tained as^ 

(S z ) = {{S~P){1 + P) 2S+1 + (l + S + P)P 2S+1 } 

x {(l + P) 25+1 -P 2S+1 }-\ (A.2) 



where P = (l/N)J2 q n ( UJ q )- The transverse two- 
spin correlation functions Cjj + are calculated from 
Eq. (|A.1[) for n = which yields 

cr + ^E«K)^ (a.3) 

The transverse correlation length £ H is calculated from 
the long-distance behavior of Eq. (|A.3j) according to 
Eq. (144")) . For comparison, the correlation length £+~ 
may be obtained from the expansion of the static spin 
susceptibility x q ~ around q — (cf. Sec. HVBh We get 




(A.4) 



The longitudinal correlation functions cannot 
be obtained by the RPA, except for the NN correla- 
tion function c[^ zz which we evaluate proceeding as 
in Ref. [l?] for S = 1/2. That is, we calculate the 
internal energy in RPA starting from the exact rep- 
resentation (f2"5|) and inserting the RPA results (|A.1|) 
and jX2|, = (l/iV)M« + - Y, q n(oJ q ) cosq x with 

M«+- = 3((5 Z ) 2 ) - (S z ) - S(S + 1), and ((S z ) 2 ) = 
S{S + 1) - (S z )(l + 2P) resulting from Eq. ©. More- 
over, we perform the decoupling C{^ zz = {S z ) {(S z ) 2 ) . 
From u, h and (S z ), the correlator C^ zz may be 

calculated. 
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